Genetic association of ovarian cancer in Asian populations
Dennis Hazelett
12/7/2017
Functional enrichment of associated variants in epigenetic data sets
Data import & munging
The main goal of this first section is to import two datasets–a foreground and background SNP list–and convert them to GRanges objects friendly to funcivar analysis. The foreground consists of Bayesian False Discovery Probability <0.5 SNPs from the study from X regions, and the background list includes all imputed SNPs of the surrounding 1MB at those regions. These lists were obtained from Johnathan Tyrer at U. Cambridge, UK.
Let’s get started with the foreground list. We need to import the list of names and genome interval locations from a file called “bfdp05.list” and transform it into a GenomicRanges object.
fg <- read.delim(file="bfdp05.list", header = TRUE) # foreground list
fg <- mutate(fg, chrom_label=paste("chr", Chromosome, sep=""))
fg$chrom_label[fg$chrom_label=="chr23"] <- "chrX"
fg.gr <- GRanges(seqnames = fg$chrom_label,
ranges = IRanges(start = fg$Position,
end = fg$Position),
strand = "*")
fg.gr$blockStarts <- start(fg.gr)
fg.gr$blockStarts <- start(fg.gr)
fg.gr$blockSizes <- end(fg.gr)-start(fg.gr)+1
## for bdfp < 0.5 SNPs "effect"=allele1, "baseline"=allele2
values(fg.gr) <- DataFrame(label=fg$SNP, linkage_group=fg$Chromosome,
allele1=fg$Effect,
allele2=fg$Baseline,
blockStarts=start(fg.gr),
blockSizes=end(fg.gr)-start(fg.gr))Next we’ll do similar for the backgroung list. The background variants consist of all the variants within 1 MB flanking on either side of the primary signal for which associations were tested. This list is kept in a file called “Selected_SNP_regions.txt”. To format it for funciVAR we need to import it and recreate it as a GenomicRanges object:
bg <- read.delim(file = "Selected_SNP_regions.txt", header = TRUE) # background list
bg <- mutate(bg, chrom_label=paste("chr", Chromosome, sep = ""))
bg$chrom_label[bg$chrom_label=="chr23"] <- "chrX"
bg.gr <- GRanges(seqnames = bg$chrom_label,
ranges = IRanges(start = bg$Position,
end = bg$Position),
strand = "*")
values(bg.gr) <- DataFrame(label=bg$SNP, linkage_group=bg$Chromosome,
allele1 = bg$A1,
allele2 = bg$A2,
blockStarts=start(bg.gr),
blockSizes=end(bg.gr)-start(bg.gr))
background.bed <- list(seqnames = seqnames(bg.gr), starts = start(bg.gr), ends = end(bg.gr), names = mcols(bg.gr)$label)
cat("browser position chr1:157878955-158078955\ntrack name='ASIAN BKG SNPS' description='Background SNPs for Ovarian Cancer study in Asian population'\n", file = "asian_oc_background.hg19.bed")
write.table(background.bed, file = "asian_oc_background.hg19.bed", sep = '\t', row.names = FALSE, col.names = FALSE, quote = FALSE, append = TRUE)Sanity check
Let’s check whether any of the foreground SNPs overlap background, i.e. are they a subset. Then we’ll assign them to a variants object.
overlaps.mat <- findOverlaps(fg.gr, bg.gr)
length(queryHits(overlaps.mat))## [1] 1313length(fg.gr)## [1] 1283length(bg.gr)## [1] 1836907unique(seqnames(bg.gr))## [1] chr1 chr2 chr3 chr4 chr6 chr7 chr8 chr9 chr10 chr12 chr13
## [12] chr14 chr16 chr18 chr20 chr21 chrX
## 17 Levels: chr1 chr2 chr3 chr4 chr6 chr7 chr8 chr9 chr10 chr12 ... chrXaoc.variants <- list(fg=fg.gr, bg=bg.gr)From this it appears that yes, foreground is a subset of background. There are a greater number of queryHits than query rows because some variants are larger than one nucleotide and can overlap more than one row in the subject.
Annotations
The annotations are mostly in-house generated epigenomics data from prior studies, rather than the StateHub treatments of Roadmap data that we typically use. I import these using the “GetSegmentations()” functionality of funciVar.
annos <- list.files("biofeatures", full.names = TRUE)
biofeatures <- GetBioFeatures(files=annos, genome="hg19") # hang on! this will take a whileEnrichment analysis
Global enrichment
First we’ll execute the standard global enrichment analysis that funciVar does, plotting enrichment values for each set of epigenetics data on and individual basis.
enrichment <- CalculateEnrichment(aoc.variants, biofeatures, feature.type="biofeatures", CI = 0.95, return.overlaps = TRUE)
#save(aoc.variants, biofeatures, file="asiandata.Rda")
global.plot <- PlotEnrichment(enrichment$enrichment, value = "difference", color.by = NULL) +
ggtitle("Global Enrichment")
#ggsave(filename = "plots/global_enrichment_plot.pdf", plot = global.plot, device = "pdf", width = 18, height = 12, units = "in")
#ggsave(filename = "plots/global_enrichment_plot.png", plot = global.plot, device = "png", dpi = 300, width = 18, height = 12, units = "in")Next we’ll create our overlap matrix for export as a supplementary table (ST7).
Locus-by-locus enrichment
Next we want to see whether individual loci are enriched for particular samples, using the same background SNPs as in the global analysis. We’ll have to define the extent of each locus and then subset the SNP lists so that we can repeat a funciVar analysis on this smaller set. The individual loci are defined in a file called “regions.txt”
locus <- with(read.delim(file="regions.txt", sep="\t"), GRanges(seqnames = seqnames, ranges = IRanges(start = start, end = end), region = region))
#define a function to create an enrichment analysis for all SNPs that overlap an interval in the locus table
l<-dim(enrichment$enrichment)[1]
w<-length(locus)
enrichment.mat <- list(probability = data.frame(matrix(NA, nrow = l, ncol = w+1)), # +1 extra column=global
difference = data.frame(matrix(NA, nrow = l, ncol = w+1)))
rownames(enrichment.mat$probability) <- rownames(enrichment$enrichment)
colnames(enrichment.mat$probability) <- c("global", as.character(locus$region))
rownames(enrichment.mat$difference) <- rownames(enrichment$enrichment)
colnames(enrichment.mat$difference) <- c("global", as.character(locus$region))
enrichment.mat$probability$global <- unlist(enrichment$enrichment$probability)
enrichment.mat$difference$global <- unlist(enrichment$enrichment$difference)
for (a in 1:w) {
print(paste(as.character(locus[a]$region)[1], "Enrichment", sep=" "))
locus.variants <- list(fg=fg.gr[subjectHits(findOverlaps(locus[a], fg.gr))],
bg=bg.gr[subjectHits(findOverlaps(locus[a], bg.gr))])
a.locus <- CalculateEnrichment(locus.variants, biofeatures, feature.type="biofeatures", CI = 0.95, return.overlaps = FALSE)
enrichment.mat$probability[rownames(a.locus),a+1] <- a.locus$probability
enrichment.mat$difference[rownames(a.locus),a+1] <- a.locus$difference
a.plot <- PlotEnrichment(a.locus, value = "difference", color.by = NULL) + ggtitle(paste(locus[a]$region[1], "enrichment", sep=" "))
ggsave(filename = paste("plots/allmarks/", locus[a]$region[1], ".pdf", sep=""), plot = a.plot, device = "pdf", width = 18, height = 12, units = "in")
ggsave(filename = paste("plots/allmarks/", locus[a]$region[1], ".png", sep=""), plot = a.plot, device = "png", width = 18, height = 12, units = "in")
}## [1] "1q23.1 Enrichment"## [1] "2p25.2-25.1 Enrichment"## Warning in SetOverlaps(variants$fg, features): No overlaps were found for
## this variant/feature combination
## [1] "2p24.2 Enrichment"## [1] "2p16.1-p15 Enrichment"## Warning in SetOverlaps(variants$fg, features): No overlaps were found for
## this variant/feature combination
## [1] "2q37.3 Enrichment"## [1] "3p25.3-2 Enrichment"## Warning in SetOverlaps(variants$fg, features): No overlaps were found for
## this variant/feature combination
## [1] "3p24.2-1 Enrichment"## Warning in SetOverlaps(variants$fg, features): No overlaps were found for
## this variant/feature combination
## [1] "3p23-22.3 Enrichment"
## [1] "4p15.33 Enrichment"
## [1] "6p25.2-1 Enrichment"
## [1] "6p21.33-31 Enrichment"
## [1] "6q26-27 Enrichment"
## [1] "7p21.3 Enrichment"## [1] "7q36.3 Enrichment"## Warning in SetOverlaps(variants$fg, features): No overlaps were found for
## this variant/feature combination
## [1] "8p12 Enrichment"## Warning in SetOverlaps(variants$fg, features): No overlaps were found for
## this variant/feature combination
## [1] "8q24.21 Enrichment"
## [1] "9q22.1-2 Enrichment"
## [1] "10q25.2 Enrichment"
## [1] "12q24.32-33 Enrichment"
## [1] "13q14.2 Enrichment"## [1] "14q31.3-32.11 Enrichment"## Warning in SetOverlaps(variants$fg, features): No overlaps were found for
## this variant/feature combination